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Abstract. Recently there has been substantial interest in least-squares finite element methods for velocity-vorticity- 
pressure formulations of the incompressible Navier-Stokes equations. The main cause for this interest is the fact that 
algorithms for the resulting discrete equations can be devised which require the soluiton of only symmetric, positive 
definite systems of algebraic equations. On the other hand, it is well-documented that methods using the vorticity as 
a primary variable often yield very poor approximations. Thus, here we study the accuracy of these methods through 
a series of computational experiments, and also comment on theoretical error estimates. It is found, despite the failure 
of standard methods for deriving error estimates, that computational evidence suggests that these methods are, at the 
least, nearly optimally accurate. Thus, in addition to the desirable matrix properties yielded by least-squares methods, 
one also obtains accurate approximations. 


1. INTRODUCTION 

The approximate solution of the Navier-Stokes equations of incompressible flow has received 
tremendous attention from engineers and mathematicians; see, e.g., [9], [10], or [11]. Among the 
more recent developments has been the use of least- squares ideas; see, e.g ., [7] for a recent survey 
of one such approach. Also, truly least-squares methods have been developed and applied; see, e.g,, 
[3]-[5], [13]-[17], [20], and [22]. 

Here, a finite element method based on a least-squares variational principle is examined for the 
approximate solution of the stationary, incompressible Navier-Stokes equations. These equations are 
cast into a first-order system of partial differential equations involving the velocity, vorticity, and 
pressure as dependent variables. In three- dimensions one has seven unknown scalar fields . However, 
the application of a least-squares principle along with, for example, a Newton linearization, results in 
symmetric , positive definite linear systems , at least in a neighborhood of the solution. The influence 
of the Reynolds number on the positive definiteness of these linear systems is felt only through the 
size of the neighborhood. Thus, if properly implemented continuation (with respect to the Reynolds 
number) methods are used, one can expect to only encounter symmetric, positive definite linear 
systems in the solution procedure. A further advantage of this method is that a single piecewise 
polynomial finite element space may be used for all test and trial functions , i. e. , one may use equal 
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order interpolation with respect to a single gnd for all dependent variables and test functions. A 
final advantage resulting from the use of a least-squares principle is that, unlike some other methods 
involving the vorticity, no artificial numerical boundary conditions for the vorticity need be devised. 

The method discussed here is similar to the ones of [3]-[5] and [13]-[17]; however, there are 
certain differences as well, especially in the order in which the least-squares, the discretization, and 
the linearizations steps are taken. Furthermore, the analyses found in some of these papers are 
incorrect, leaving open the question of the accuracy of approximations. 

In §2, we define the least-squares finite element method. In §3, we discuss some practical and 
theoretical issues connected with the method described in §2. Then, in §4, we give the results of 
a computational study of the accuracy of the algorithm. Finally, in §5 we give some concluding 

remarks. 


2. THE LEAST-SQUARES FINITE ELEMENT METHOD 

2.1 - The velocity- vorticity-pressure equations 

Let the bounded set ft C R 3 denote the flow domain and let T denote its boundary. The 
dimensionless equations governing the steady incompressible flow of a viscous fluid may be written 
in the form 

(2.1) divu = 0 in ft , 


( 2 . 2 ) 


curl u — u> = 0 in ft , 


(2.3) 

and 


i/curlw + u> x u + gradp = f in ft, 


(2.4) divu> = 0 in ft, 

where u, p, and a> denote the velocity, pressure, and vorticity fields, respectively, v the inverse of the 
Reynolds number, and f a given body force. The first-order system of partial differential equations 
(2.1)-(2.4) constitutes the velocity-vorticity-pressure form of the equations of steady, incompressible 
flow. Note that the terminology “pressure” is a little misleading since the variable p is actually the 
total head, i.e., p = p+ ^|u| 2 , where p denotes the pressure. 

In view of (2.2), it seeems that (2.4) is redundant. However, it is crucial to the stability of the 
least-squares algorithm to explicitly require (2.4); see, e.g., [3] and [15]. 

The system (2.1)-(2.4) should be supplemented with boundary conditions. Here we examine 
two such boundary conditions. The first is to impose the velocity on the boundary, i.e., 

(2.5) u = U x on T , 
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where Ui denotes a given function defined along T. Note that (2.5) implies that the normal com- 
ponent of the vorticity is also known, i.e., 

(2.6) a; • n = n * curlUi onT, 

where n denotes the unit outer normal to fl. To see this, one merely needs to observe that n * curl u 
involves only tangential derivatives of the tangential components of u and these may be deduced 
from (2.5). To these one must add a condition to fix the pressure; we choose to fix the pressure at 
a single point xo in Cl, i.e. y 

(2.7) p(xo) - Po , 
where Po is a given number. 

The second combination of boundary conditions we consider is the pressure and the normal 
component of the velocity, i.e., 

(2.8) p = P on T 
and 


(2.9) u n = J7 2 onT, 

where P and C/ 2 denote given functions defined along T. The boundary conditions (2.8)-(2.9) are not 
all that useful as boundary conditions for the Navier-Stokes equations; we consider them here for 
two reasons. First, a complete, rigorous analysis of least-squares finite element approximations of 
(2.1)-(2.4) and (2.8)-(2.9) can be given using standard techniques; this is not the case for (2.5)-(2.7). 
Furthermore (2.1)-(2.4) and (2.8)-(2.9) can be shown to be related to second order elliptic partial 
differential equations; we will discuss these issues in more detail below. 


2.2 - The least-squares principle 

One can use (2.1)- (2.4) to define the least-squares functional 


( 2 . 10 ) 


J7(u,u>,p) = J ^jcurlu- u >\ 2 -f |divu| 2 + |diva;| 2 

+ |z^curl u) + uj x u 4- gradp - f | 2 ^ <Kl . 


The least-squares principle then requires the minimization of this functional over appropriate func- 
tion spaces. Then standard techniques from the calculus of variations may be used to deduce that 
minimzers (u,p,o>) of J necessarily satisfty 


( 2 . 11 ) 


/ (curl u — a>) • curl v -f div u div v 

+(^curlu; + gradp + u) x u — f) • (o> x v)j d£l = 0 , 
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(2.12) 

J ^divu? div£ — (curlu - w) ■ f 

n A 



+ (i/curlu; + gradp + u> x u - f) * (1 / curlC + C x U )J 

and 



(2.13) 


/ (t/curlo; + gradp + a>xu-f) - gradgdfi = 0 . 

Jci 


In (2. 11)- (2. 13) the test functions (v, q, £) are required to belong to suitable function spaces; we do 
not go into detail here since we are primarily interested in finite element discretizations of these 
equations. 

Clearly any solution (u,w,p) of, say (2.1)-(2.7), satisfies (2.11)-(2.13). 


2.3 - The two-dimensional case 

For planar flows we have that u = (u u in, 0) T and u u u 2 , and p are functions of x\ and x 2 
only. Then, we have that w = (0, 0, w)^ where w = dinj dx\ — du\ j dx 2 . In this case, the system 
(2.1)- (2.4) simplifies to 


(2.14) 


du\ du 2 
dxi + dx 2 


in Cl , 


(2.15) 


du 2 

dx\ 


- u = 0 in Cl , 
dx 2 


(2.16) 


du> dp , 

!/-— + u 2 w = f 1 

OX 2 dx I 


in Cl , 


and 

(2.17) 



dp 

+ a - +UiW = 

0X2 


h 


in Cl, 


where f = (/i, h . °) T - The boundary conditions (2.5)-(2.7) reduce to just (2.5) and (2.7) and 
the boundary conditions (2.8)-(2.9) remain unchanged. The functional (2.10) and the necessary 
conditions (2.11)-(2.13) also simplify in the obvious manner for two-dimensional problems. 


2.4 - Finite element methods 

Starting with the weak formulation (2.11)-(2.13), a conforming finite element method can be 
defined in a completely standard manner. We choose a finite element space S h parametrized by 
h. For example, for a given positive integer r, S h could consist of continuous (over Cl) piecewise 
polynomials of degree less than or equal to r with respect to a subdivision of Cl into finite elements. 
In this case the parameter h may be related to the size of the grid. We then define the spaces 
S h = {v\vi€S h ,i = 1,2,3}, 

Vo = {v€S h |v = 0onr}, 
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ZS = {C€S 7 , |C-n = Oonr} ) 


and 


Qo = { q£S h | g(xo) = 0 }. 


For the boundary conditions (2.5)-(2.7), the discrete problem is defined as follows: seek u h £ S h , 
w* £ S h , and p h £ S h such that u h = Uf and u h • n = W h on T, p fc (x o) = Po, 


(2.18) 


J J(curl u h - w h ) ■ curl v h + div u h div v h 

+(«/curlw h + gradp h + w h x - f) ■ (ut h x v h )j d£l = 0 V V* € Vq , 


(2.19) 


J [divu>* div^ - (curl u h - u h ) • £* 

+ (i/curlw 7 * + gradp h + u h x u h — f) • (i/curl^ 1 + C* x u h )j <Kl = 0 V £* £ Zq , 


and 

(2.20) f (i/curlw h + gradp** + u> h x u h — f) ■ grad q h d£l = 0 V q h £ Qq 

Jo 

are satisfied. Here U* and W h are approximations of the data Uj and curlUi • n, respectively. For 
example, we could define the former pair to be boundary interpolants of the latter pair. 

Note that all of the discrete variables, i.e., q h and the components of u h and U3 h , are approxi- 
mated by the same degree continuous piecewise polynomials defined with respect to a single grid. 

The discrete problem for the boundary conditions (2.8)-(2.9) is also given by (2.18)-(2.20) except 
that now we have that u h • n = t/j and p h = P h on T and the test functions spaces in (2.18)-(2.20) 
must be suitably redefined to account for the boundary conditions v A • n = 0 and ^ = 0 on T. 
Again, Uff and P h are approximations to U2 and P, respectively. 


3. PRACTICAL AND THEORETICAL ISSUES 
3.1 - Newton’s method 

The discrete equations (2.18)-(2.20) are a nonlinear system of algebraic equations that must 
be solved in an iterative manner. There are many methods that one might use for such a purpose; 
here we only consider Newton’s method. However, we do remark that if a quasi-Newton method is 
used, it should be chosen so that it preserves symmetry and positive definiteness of the approximate 
Hessian matrices. 

Newton’s method for the solution of (2.18)-(2.20) is defined as follows. Given initial guesses u^°l, 
cj(°\ and p(°l for u h , u> h , and p h , respectively, the sequence of Newton iterates { p^ }k>o 
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is generated recursively by solving, for k = 1,2, . . the system 


(3.1) 



. curlv fc + divu^ divv h 


+ (i/curlo; (fc) + gradp (fe) + w (fe) x u (k_1) + u> (fc_1) x u (fc) ) • (w (fc x) x v h ) 
+ (j/curla^-^ + gradp^ + J k ~ 1) x u< fc - x > - f) • (w (fc) x v h )] dfi 


J (t/curlw (fc_1) + gradp (fc_1) + 2w (fc x) x u (fc x) ) • (u> (fc x) x v h ) dZl V v h € V£ , 


(3.2) 


J j^divt*/*^ div^ 7 * — (curlu^ 




+ (j/curl w (fe) + gradp (k) + w (fc) x u (fc_1) + u (k x) x u (fe) ) • (i/curlC 7 *) 

+ (t/curlw^ + gradpW +»« x W*" 0 x u' fc ’) • «* x u^' 1 *) 

+ (i/curlw (fc-1) + gradp^- 1 ) + w* 1 - 1 * x u^ 1 * - f) • (C* x u«)] dfi 


= J^v(f + u>( k ^ x u^ fc J )) • curlC 7 * 

+ (r/curlw^-^ + gradp< fc - 1 ) + 2w< fc - 1) x *<*-») - (<* x n^^)] dft V C* € Z* , 


and 


(3.3) 


f (t/curla>( fc ) + gradp^ +u^ x u* fc ^ +<*/* ^ 

= f (i + u) ( - k ~ 1 ^ x u^ fe-1 ^) -grad^tKl 

Jn 


Ja 


x u«) • gradg k dfl 


The system of linear algebraic equations (3.1)-(3.3) that is used to determine the fc-th Newton 
iterate from the (fc — l)-st looks rather formidable. However, it also has some very good features. 
First, it is easy to see that this system is symmetric. Moreover, in a neighborhood of a minimizer, 
the Hessian matrix for the functional of (2.10) is necessarily positive definite; but this Hessian matrix 
is exactly the coefficient matrix of (3.1)-(3.3). Thus, in a neighborhood of a solution of (2.18)-(2.20), 
the system (3.1)-(3.3) is symmetric and positive definite. This feature is independent of the value of 
v, i.e., of the value of the Reynolds number. These observations, along with the guaranteed local 
and quadratic convergence of Newton’s method, are potentially very advantageous. 


3.2 - Continuation methods 

At this point one may well ask what goes wrong with the method as the Reynolds number 
increases? Surely one difficulty is that the attraction ball for Newton’s method, or for any other 
iterative method for solving the nonlinear equations, decreases in size with increasing Reynolds 
number. This is known to be true for other discretizations of the Navier-Stokes equations. A 
related observation is that the positive definiteness of the Hessian matrix is guaranteed only in a 
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neighborhood of a minimizer; again the size of this neighborhood surely decreases with increasing 
Reynolds number. As a result, for an arbitrary initial guess we may have that Newton’s method 
doesn’t converge and/or that the coefficient matrix in (3.1)-(3.3) is not positive definite. The former 
is, of course, unacceptable, while the latter occurance would preclude the use of simple iterative 
methods for solving the linear systems (3.1)- (3.3) that define the Newton iterates. 

In order to determine an initial guess that is within the attraction ball of Newton’s method, 
and also is such that the coefficient matrix in (3.1)-(3.3) with fc = 1 is positive definite, one can use 
continuation or homotopy methods, among others. Here we describe a simple continuation method; 
see, e.#., [18], [19], or [21]. Let us symbolically express the system (2.18)-(2.20) in the form 


F(u h ,u h ,p h ;Rej =0. 


Here Re = 1 jv denotes the target Reynolds number, i.e., the value of the Reynolds number at which 
we want a solution. Now, suppose we have a sequence of increasing Reynolds numbers {fle m }^ =1 
with Re\f = Re . We denote the solution of (2.18)-(2.20) for v m = 1/Rem by (u m ,w m ,p m ). For 
any value of m, we obtain (%, cj m ,p m ) by Newton’s method, *.e., we solve the sequence of linear 
systems 

■ (<■£?. ■.«>,?£>) - (itf- 1 ’, jfi- 1 ’)) 

<3 ' 4> = -f (u<;- 1 >,<.£- 1) .Pm" 1) ; fle m) for k = 1,2 

Here, F' denotes the Jacobian of F with respect to (u ft , u^,p h ). We need to specify the initial guesses 
(um\wS\p^) to be used to start, for each m, the iteration with respect to A: in (3.4). Assume 
that Re i is sufficiently small so that the iteration (3.4) converges if (u^ 0) ,i«/f > ,Pj 0) ) is chosen to 
be the solution of a discrete linear Stokes problem. For example, we could choose Re i = 1 and 
take for (uf^wf^p^) the solution of (2. 18)- (2.20) with all cross product terms deleted and with 
v = 1. The latter problem involves a symmetric, positive definite linear algebraic system. The 
remaining initial guesses {(u^ ) ,a>^ ) ,p^ ) )}^f =2 are determined by “continuing along the tangent”, 
i.e., by solving the linear algebraic system 

Rem-l) ' ((“S^m >Pm ) " K-l.“>»-l.Pm-l)) 

= {Ren Re m — l)FRe^(Um — lj — 1 1 Pm — 1 1 R&m — l) ■ 

Here, (um-^Wm-i.pm-i) denotes the converged iterates determined from (3.4) at the (m — l)-st 
stage and Fut denotes the Frechet derivative of F with respect to the parameter Re. Note that the 
coefficient matrix of the linear system (3.5) may be chosen to be the same as the coefficient matrix 
for the last iteration of (3.4) at the (m - l)-st stage. 

The combined Newton-continuation method is now completely defined. If {Re m - Re m -i) is 
sufficiently small, the use of (3.5) should yield initial guesses that are within the attraction ball of 
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Newton’s method and such that the coefficient matrices in (3.4) with k = 1 are positive definite. 
(They are always symmetric, of course.) In fact, since Newton’s method is guaranteed to be locally 
convergent, i. e., its attraction ball is nontrivial, and since the neighboorhood of a minimizer for which 
the Hessian matrix is positive definite is also nontrivial, by choosing (Re m —Rem- 1 ) sufficiently small, 
one can guarantee that the combined Newton-continuation method should only deal with symmetric, 
positive definite matrices. 

The method can be made self-correcting. For example, suppose that the linear systems in (3.4) 
are solved by an iterative method, e.g., the conjugate gradient method, that works, or works well, 
only for symmetric positive definite matrices. Then, if we have chosen an increment (R&m — R^m-i) 
that is too large, then either the Newton iteration or the conjugate gradient iteration will fail. In 
either case, one can restart the m-th stage by choosing a smaller value of Rem in (3.4) and (3.5). 

We note that often the even simpler “continuation along a constant” method 


(3.6) •J’S?) = (u m -l,«m-l,Pm-l) 

for generating initial guesses has been found to work well in viscous flow calculations; see [12]. The 
disadvantage of (3.6) it that it breaks down in the vicinity of bifurcation or turning points, while 
(3.5) can be ammended so that it can handle such angular points; see [18], [19], and [21]. 

We also note that there are modifications possible to Newton’s method that circumvent the 
difficulty of not having positive definite Hessian matrices. For example, one ample such modification 
is to add to the diagonal entries of the coefficient matrix in (3.1)-(3.3) a multiple of the magnitude of 
the residual of the previous Newton iterate. In the notation of (3.4), we would replace the Jacobian 
matrix F' on the left-hand side with 


(3.7) 


where || • || denotes the Euclidean length. By choosing the constant 7 sufficietnly large we can make 
sure that the matrix represented by (3.7) is positive definite. However, as we approach the solution 
of the problem, the term multiplying 7 approaches zero so that we recover the local convergence 
properties of Newton’s method in the neighborhood of the solution. 


3.3 - Enhancing mass conservation 

In many instances one is especially interested in conserving mass, t.e., satisfying the continuity 
equation (2.1) as well as possible. The method discussed here is not exactly mass conserving, i.e., 
divu h ^ 0 exactly. In fact, one can easily show that at best divu h « Ch r whenever the finite 
element space, restricted to any element, contains all polynomials of degree less than or equal to r. 
We can reduce the size of the constant C by introducing a weight into the functional (2.10). Indeed, 
if we replace the |divu| 2 term in (2.10) by a|divu| 2 , where a > 0 is a constant, then the constant 
C can be shown to be proportional to 1/y/a. Thus, by choosing a large value of a, one can make 
divu h small. This is another potential advantage of the least-squares method. However, one must 
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keep in mind that the larger the value of a, the worse one satisfies the momentum equation relative 
to the continuity equation. 

3.4 - Theoretical observations concerning accuracy 

Error estimates for the least-squares finite element approximations are derived through the 
following process. First, using the theory of [2] (see also [6] and [9]), one can show that the error 
estimates for the nonlinear Navier-Stokes equations are essentially the same as those for the linear 
Stokes problem, at least away from singular points. Now for the Stokes problem 


(3.8) 

divu =0 in fl, 

(3.9) 

curl u — u = 0 in , 

(3.10) 

z^curlo; + gradp = f in Q , 

and 


(3.11) 

diva; =0 in 


standard techniques for estimating the error of least-squares finite element approximations require 
that this system be elliptic and that the boundary conditions satisfy the complementing condition 
of [1]; see, e.g., [4], [5], and [23]. Moreover, in order to use finite element functions that are merely 
continuous across element boundaries, one must be able to bound the L 2 -norm of the derivatives 
of solutions of (3.8)-(3.11) in terms of the data of the problem. For the boundary conditions (2.5)- 
(2.7) this program cannot be carried out, i.e., the system (3.8)-(3.11) and (2.5)-(2.7) does not 
satisfy the complementing conditions of [1] when one requires that the derivatives of u, o>, and p be 
merely square integrable. Thus, standard error estimation techniques for least-squares finite element 
discretizations of (2.1)-(2.7) cannot be used. 

On the other hand, one can easily show that the system (3.8)-(3.11) and (2.8)-(2.9) does satisfy 
all the requisite conditions of the theory of [1] and therefore one may conclude that in that case one 
may obtain optimal error estimates. For example, if continuous piecewise polynomials of degree < r 
are used with respect to a triangulation of Q, then one can ultimately conclude that for sufficiently 
smooth solutions 

(3.12) |u- u h |i + |o> -o^li + |p-p'Mi < Ch r 
and, under mild restrictions on the domain, 

(3.13) || u - u h || 0 + Ik -a^Ho + \\P -/Ho < Ch r+1 , 


9 


where h denotes a measure of the grid size and where || • ||o and | • |i denote the L 2 (Q)-norm and the 
L 2 (ft)-norm of the first derivatives, respectively. 

Nothing we have said necessarily means that least-squares finite element approximations of 
(2.1)-(2.7) do not satisfy (3.12)-(3.13) as well. All we have concluded so far is that we cannot, 
using standard estimation techniques, derive these estimates. Indeed, this is the motivation for the 
computational study of §4. 

3.5 - Connection of (2.8)-(2.9) with second-order elliptic problems 
Consider the problem 

(3.14) &P = 9 infi, 
where A denotes the Laplacian operator, and 

(3.15) P = P on r . 

Now, let 

(3.16) u = gradp in ft. 

Then we have that 

(3.17) divu = p in ft. 

Thus (3.16)-(3.17) is a first-order system equivalent to (3.14). However, it can be shown that this 
first-order system is not elliptic and that in general least-squares finite element methods for (3.16)- 

(3.17) are not stable; see [8]. However, if one considers the system (3.17), 

(3.18) curl u> + gradp = u, curl u = 0 and div<*> = 0 in ft, 

then it can be shown that this sytem is elliptic. Clearly, if p is a solution of (3.17)-(3.18) and (3.15), 
then p is also a solution of (3.14)-(3.15). Note that the principal part of (3.17)-(3.18), i.e., the 
differentiated terms, are identical to those of (2.1)-(2.4), and that the boundary condition (3.15) is 
the same as (2.8). We can then infer the optimal accuracy of least-squares finite element methods 
for appropriately defined first-order systems arising from second-order ellliptic equations. 

4. A COMPUTATIONAL STUDY OF ACCURACY 

The accuracy of least-squares finite element approximations for the Navier-Stokes equations are 
essentially the same as that for the linear Stokes equations. Moreover, the accuracy is independent 
of dimension. Thus, here we restrict attention to the Stokes equations in two-dimensions. We also 
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take for our domain the unit square £1 = { 0 < x < 1 , 0 < y < 1 }. Thus, we consider the generalized 
Stokes equations 


(4.1) 


du dv 
dx + dy Sl 


in 


(4.2) 


dv 

dx 


— --w = g 2 in S2, 
dy 


(4.3) 


0p , 

dy dx 11 


in £2, 


and 

(4.4) 


du> dp 
dx dy 


= h 


in ft , 


where gi, gi, /i, and fa are given functions. We consider the two sets of boundary conditions 


(BC1) 


u — U and v = V on T and p(0, 0) = Po 


and 


(BC2) 


p = P and uni + v n 2 — W on V , 


where U, V, P, and W are given functions defined on r, P 0 is a given number, and ni,n 2 denote 
the components of the unit outer normal. We will define the various data functions by choosing an 
exact solution (u, v t u>,p) and then substituting into the above equations. 

We will measure the differences 


(L 2 error) 


IK-^llo 



and 

( H 1 error) 








where f could be any of u , v, a;, or p. 

Throughout we use piecewise quadratic finite element spaces; thus, for sufficiently smooth solu- 
tions and for the domain Q = {0<x<l,0<y<l}, we expect that for the boundary conditions 
(BC2) we have that 

(4.5) U-Z h \\o = 0(h 3 ) and \£ - £ h \x = 0(h 2 ), 
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where again £ could be any of u, v, t<;, or p. This is confirmed by the computations that follow. 
Of special interest to us here is the computationally determined rates of converge for the boundary 
conditions (BCl). 

The first example we present has the smooth exact solution 


(4.5) 


u = — cos ttx sin 7ty + 1 — y 3 , v = sin 7rx cos ny 4 * 1 — x 3 > 

dv du , . ,2 

v s — , and p — smycosx + xy . 

OX cfy 


Figures 1 and 2 give log-log plots of the L 7 and H l errors vs. the number of grid intervals in 
each direction for a uniform grid spacing. (Similar results have been obtained for nonumform grid 
spacings.) From these figures one can confirm that for the boundary conditions (BC2) one does 
indeed obtain the best approximation convergence rates of (4.5). (As is usually the case, computed 
L? rates are usually less reliable than their H * counterparts.) In fact, the asymptotic slopes of 
each of the dotted curves of Figure 1 are approximately 3 while those of Figure 2 are approximately 
2. Surprisingly, the slopes of the solid curves are asymptotically the same as those of the dotted 
curves; thus, for the boundary conditions (BCl) we seemingly obtain the best approximation rates 
of convergence of (4.5). Experiments with other smooth solutions yield similar results. 


L2 Error for U; BC1/BC2 L2 Error for V; BC1/BC2 






Fig. 1. I ? errors vs. number of grid intervals in each direction for example (4.5). 
Solid line if for BCl; dotted curve is for BC2. 
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Hi Error for U; BC1/BC2 


HI Error for V; BC1/BC2 






Fig. 2. /f 1 errors vs. number of grid intervals in each direction for example (4.5). 

Solid line if for BC1; dotted curve is for BC2. 

The next set of examples uses the exact solution 

(4.6) u = v = lj = p= ((x - a) 2 + (y- b) 2 ) s/2 

We use a = b = 0.1234. The exponent s may be chosen to adjust the smoothness of the exact 
solution, and thus the rates of convergence of the best approximations out of the finite element 
space used. Thus, except for some exceptional integer values, we have that 

(4.7) ||C - i h ||o = 0(h t+l ) and |f - | h |i = 0(h*) , where t = min(s, 2) , 

and where again £ denotes any of u, v, u>, or p, and £ h denotes the corresponding best approximation. 

Figures 3-10 give log-log plots of L 2 and H 1 errors vs. the number of grid intervals in each 
direction for a uniform grid spacing and for different values of s. According to (4.7), the best 
approximation for the cases s = 2.5 and s = 2.001 have L 2 and errors of 0(h , 3 ) and 0{h 2 ), 
respectively. Figures 3-6 show that the least-squares finite element solutions for both boundary 
conditions (BCl) and (BC2) also seem to converge at approximately this rate for these values of s. 
Similar conclusions can be drawn from Figures 7 and 8 for which the best approximations and both 
finite element solutions have L 2 and H 1 errors close to 0(h 2 5 ) and 0(/i 15 ), respectively. Likewise, 
from Figures 9 and 10 we have that finite element solutions have L 2 and H l errors close to 0(h 2 ) 
and 0(/i 1 ), respectively; these rates are again those of the best approximation for s = 1.001. 
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Fig. 4. H l errors vs. number of grid intervals in each direction for example (4.6) with s = 2.5 
Solid line if for BC1; dotted line is for BC2. 
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L2 Error for U; BC1/BC2 



L2 Error for W; BC1/BC2 



L2 Error for V; BC1/BC2 


0.00005 

0. 0000JL 
5. 10 

-6 

1. 10-7 
5. 10 



23 57 10. 15. 


L2 Error for P; BC1/BC2 

0.001 


0.0001 



0.00001 


-6 


1. 10 



23 5 7 10. 15. 


Fig. 5. I? errors vs. number of grid intervals in each direction for example (4.6) with s = 2.01 
Solid line if for BC1; dotted line is for BC2. 


HI Error for U; BC1/BC2 



HI Error for V; BC1/BC2 



HI Error for W; BC1/BC2 


HI Error for P; BC1/BC2 




Fig. 6. H l errors vs. number of grid intervals in each direction for example (4.6) with s = 2.01 
Solid line if for BC1; dotted line is for BC2. 











Fig. 8. H 1 errors vs. number of grid intervals in each direction for example (4.6) with s — 1.5 
Solid line if for BC1; dotted line is for BC2. 
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The above conclusions drawn from Figures 1-10 can be quantized by computing the slope of a 
least-squares straight line fit to the various curves in the figure. The results are given in the following 
table. The first column of the table gives the exact solution used to obtain the remaining columns. 
The second column gives the error type, i.e., either H 1 or L 2 , to which the entries of the remaining 
columns correspond. The third column gives the rate of convergence of the best approximations 
to the exact solutions out of the quadratic finite element space. Specifically, that column gives the 
value of t for the H 1 error and t + 1 for the L 2 error as determined from (4.7) The fourth column 
gives the exponent j3 in the error formula 

\\u-u h \\=0(h^, 

where || • || denotes either the L 2 norm or the H l semi-norm, u the exact x-component of the velocity, 
and u h its least-squares finite element approximation. The exponent /? is determined computationally 
from the data that was used to produce the plots for the errors in the approximation to u in Figures 
1-10. Specifically, 8 was determined as the slope of the best least-squares straight line fit of the 
curves in those figures. The remaining three columns give similar data for v, u, and p. 


Exact solution 

BC type 

Error type 


Rates of convergence 





b.a. 

u 

V 

U) 

P 

(4.5) 

BCl 

H 1 

2 

1.97 

1.96 

2.01 

2.09 


L 2 

3 

3.54 

3.61 

3.44 

3.21 


BC2 

H 1 

2 

1.98 

1.98 

1.98 

2.00 



L 2 

3 

3.10 

3.11 

3.04 

3.00 

(4.6) s = 2.5 

BCl 

H 1 

2 

1.93 

1.93 

1.88 

1.89 


l? 

3 

3.00 

2.98 

2.88 

2.92 


BC2 

H 1 

2 

1.93 

1.93 

1.90 

1.94 



L 2 

3 

2.96 

2.97 

2.89 

2.94 

(4.6) s = 2.01 

BCl 

H 1 

2 

1.98 

1.98 

1.86 

1.85 


L 2 

3 

2.82 

2.81 

2.70 

2.98 


BC2 

H 1 

2 

1.98 

1.98 

1.84 

1.76 



L 2 

3 

2.93 

2.95 

2.71 

2.87 

(4.6) s = 1.5 

BCl 

H 1 

1.5 

1.41 

1.41 

1.82 

1.82 


I? 

2.5 

2.39 

2.42 

2.52 

2.36 


BC2 

H l 

1.5 

1.39 

1.39 

1.43 

1.41 



L 2 

2.5 

2.39 

2.42 

2.52 

2.36 

(4.6) s = 1.001 

BCl 

H 1 

1.001 

0.99 

0.99 

1.34 

1.37 


L 2 

2.001 

1.87 

1.87 

1.98 

1.81 


BC2 

H 1 

1.001 

0.98 

0.98 

0.97 

0.97 



L 2 

2.001 

1.92 

1.94 

2.06 

1.89 


X&ble 1 . Rates of convergence of the and L 2 errors in the best approximation (b.a.) 
and in the least-squares finite element solution (u, V, u>, p ). 
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Note that the behavior of the L 2 errors is much more erratic than that of the H x errors; this is 
usually the case. However, we see that the errors for (BCl) are not all that different from those for 
(BC2) and we may conclude that the errors in the former case are, at the least, nearly optimal. 

5. CONCLUDING REMARKS 

We have studied the accuracy of a least-squares finite element method for the Navier-Stokes 
equations based on a velocity-vorticity-pressure formulation. In particular, we have focused on 
velocity boundary conditions. In this case standard techniques for error estimation fail so that a 
computational study of accuracy is called for. In spite of this failure, the computational experiments 
reported on here indicate that the accuracy is, at the least, nearly optimal, i.e., that the rate of 
convergence as the grid size is reduced is seemingly the same as that for the best approximations 
out of the finite element spaces. When the least- squares finite element discretization algorithm is 
coupled with a Newton linearization with continuation, desirable discrete algebraic properties result. 
Thus, the overall method seems to provide a good combination of accuracy and efficiency. 

There remains substantial issues to study connected with the least-squares finite element method 
for incompressible flows. These include practical implementation issues such as the use of iterative 
linear solvers and theoretical issues such as the derivation of rigorous error estimates. We will address 
these issues in a forthcoming paper. 
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